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Abstract. Let R = C[x] be a polynomial ring with complex 
coefficients and Dx = C(x,d) be the Weyl algebra. Describing 
the localization Rf = R[f ] for nonzero / E R as a Z?x-modulc 
amounts to computing the annihilator A = Ami(/ a ) C Dx of the 
cyclic generator f a for a suitable negative integer a. We construct 
an iterative algorithm that uses truncated annihilators to build A 
for planar curves. 

Introduction 

Let f £ R := C[xi, . . . ,x n ] be a nonzero polynomial (n being a 
positive integer) and Rf the localization of the polynomial ring R with 
respect to /. Elements in Rf are rational functions 4% with g G R 
and k G N. If / is not a constant, the -R-module Rf is not finitely 
generated. One fundamental result by J. Bernstein states that Rf is 
finitely generated when considered as a left module over the complex 
Weyl algebra A n := A n (C) of order n (see the needed definitions and 
precise statements in the next Section). This is a fundamental result 
in D-module theory (the algebraic theory of systems of linear partial 
differential equations) which can be considered as part of Singularity 
Theory. 

Bernstein's result states even more: the module Rf is cyclic over A n 
and in fact it is generated by a rational function of the form f a for some 
negative integer a. That means that Rf can be described as a quotient 
Arvn{f a ) wnere Ann(f a ) := AnriA n (f a ) is the left ideal in A n formed of 
the hnear differential operators annihilating the rational function f a . 
By [H] (see also [29] and [23]) one can take a = —n + 1. 

One main problem in algorithmic D-module theory is to compute 
a finite system of generators of the annihilating ideal Ann(f a ). There 
are several algorithms solving this problem (see e.g. [TH], [21], [20]). 
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which use Grobner bases and elimination theory in the Weyl algebra 
A n . In the worst case, computing Grobner bases has a doubly exponen- 
tial complexity in both (commutative) polynomial rings and the Weyl 
algebra (see [1] and p3]); however, in practice, one should expect much 
longer running times for the latter on the input of the same size. 

As the Weyl algebra A n is a Noetherian ring, one can associate to 
the couple (/, a) the smallest integer n(f a ) such that Ann(f a ) is gener- 
ated by operators of order less than or equal to n(f a )- This numerical 
invariant plays a relevant role in the so called Logarithmic Comparison 
Problem with respect to the hypersurface defined in the complex affine 
space C™ by the polynomial equation / = (see e.g. [TO], [27]). 

In this paper we describe a new algorithm computing Ann(/ _1 ) for 
any reduced complex polynomial / = f(x, y) in two variables. This 
algorithm uses Grobner bases techniques in both polynomial rings and 
the Weyl algebra A2, but avoids elimination theory in the latter: our 
experiments show that the bottleneck of the algorithm is a syzygy 
module computation over the former. 

We first compute iteratively, for d — 1,2, ... , truncated annihilating 
ideals Ann^(f~ 1 ) generated by linear differential operators annihilat- 
ing f~ l and of order less than or equal to d. The ideal Ann^ d \f~ 1 ) can 
be computed by using polynomial Grobner basis. We then compute 
equations for the characteristic cycle of the ideal Ann^(/ _1 ) and com- 
pare it with the one of Ann(f~ 1 ) which only depends on the multiplicity 
of the plane curve / = at its singular points; for the beforementioned 
comparison we only need to localize at each of these singular points. 

The paper is organized as follows. In the first Section we survey 
basic results on the Weyl algebra, holonomic .D-modules, 6-functions 
(or Bernstein-Sato polynomials) and the annihilating ideal of some ra- 
tional functions on C n . In Section 2 we describe the new iterative 
algorithm, we give a stopping criterion and prove its correctness. We 
also treat some examples, in particular the so called family of Reiffen's 
curves f Pyq . Our implementation of the algorithm in Macaulay2 works 
on larger examples of Reiffen's curve family than the known general al- 
gorithm and our experiment demonstrates an interesting (very simple!) 
dependence of the order of generation of Ann^(/^) on p. In the last 
Section we conclude the discussion and propose some open questions 
on the subject. 

I. Preliminaries 

1.1. Weyl algebra. Let n > be an integer and X = C n be the n- 
dimensional complex affine space. Define the n-th Weyl algebra as the 
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associative algebra 

D x = C(x, d) = C{xi, ...,x n ,di,...,d n ) 

where [<9j,£j] = d{Xi — X{d{ = 1 and all other pairs of generators com- 
mute. The Weyl algebra D x is isomorphic to the algebra of linear dif- 
ferential operators with coefficients in the polynomial ring R = C[x] = 
C[xi, . . . , x n ]. Every element in D x has a unique normal form 

Q= Qapx a d^, 

where finitely many of q a p G C are nonzero and where x a d stands for 
x" 1 ■ ■ ■ x^di 1 ■ ■ ■ d^ n . We denote the n-th Weyl algebra A n if we need 
to emphasize the dimension of X; throughout this paper Dx = A 2 is 
used the most. 

The ring Dx is simple: there are only trivial two-sided ideals. All 
.Dx-ideals and Dx-modules considered in this article are left ideals and 
modules, respectively. 

Examples of .Dx-modules include functional spaces, e.g., polynomial 
functions C[x] and smooth functions C°°(X), and power series rings 
C[[x}}. 

Another example is the localization of the polynomial ring C[x, f x ] 
where / is a nonzero polynomial with the natural action defined as 
follows: 

Xi ■ gf~ 3 = Xigf~\ 

\dXi dxij 
for 1 < i < n, g G C[x], and j G N. 

1.2. Grobner bases. It is possible to compute in the Weyl algebra 
Dx, since it is Grobner friendly; see e.g. [El [9], [21] and [13] . 

Grobner bases can be computed with respect to any ^-compatible 
monomial order, where w = (w x , wq) G K 2n satisfies w x + wq > com- 
ponentwise. This condition, in particular, guarantees that the filtration 
{F? } iez of D x , where 

F™ = Span {x a d 13 : a, f3 G N n , w (a, /3) < i) , 

is preserved under taking the commutator. Therefore, the associated 
graded algebra gr w (Dx) is well defined. Note that 

• if w x + w 9 = then gr w (D x ) = D x ; 
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• if w x + wq > componentwise then 



gr w (D x ) = C[x,^] = C[x 1 




a polynomial ring in 2n variables. 

The w-order of a nonzero operator Q = Yl a em n Qapx a d^ G Dx, 
denoted ord w (Q), is defined as the maximum of w ■ (a, (3) for q a/3 ^ 0. 
The initial part of Q with respect to w is 



in w (Q) = Yl <laf3X a ^ e gr w (D x ). 



w-{a,f})=ord w (Q) 



We simply denote in w (0) = 0. 

Definition 1. Let I be an ideal in Dx- The characteristic ideal 



with w = (0, e), i.e.: w(xi) = and w(di) = 1 for i = 1, . . . ,n defines 
the characteristic variety 



Here ia w (I) stands for the ideal of gr w (D x ) generated by {m w (Q), \ Q G 
/}, the set of the leading forms of the elements of I with respect to the 
weight w. The set V( J), for an ideal J in C[x, £], is the variety in C 2n 
defined by the ideal J. 

If M is a nonzero finitely generated _Dx~module, one can define, 
following P], its characteristic variety Char(M) C C 2n by using any 
good filtration on M. The set Char(M) is an affine algebraic subset of 
C 2n . In particular one can consider its Krull dimension dim(Char(M)). 

To each irreducible component C of Char(M) one can associate its 
multiplicity mc which is a positive integer number. For the definition 
of the Krull dimension of an affine algebraic set and the multiplicity of 
its irreducible components one can consult [I]. 

Theorem 1 (Bernstein's inequality; [2]). Let M be a nonzero, finitely 
generated Dx-module. Then n < dim(Char(M)) < 2n. In particular, 
if I £ D x is a left ideal, then n < dim(Char(Dx//)) < 2n. 

There is a number of software implementations of D-module algo- 
rithms: kan/sml [26], risa/asir [18], dmod. lib library [16] in Singu- 
lar [TTJ, and the Dmodules package [17] in Macaulay2 [12]. The last 
package hosts our implementations of the algorithms in this article. 



m w (I)CgT w (D x ) = C[x,$], 



Char(L> x /7) = V(gr ro (J)) C C 
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1.3. Holonomic D-modules. A finitely generated -D^-module M is 
called holonomic if either M = (0) or the dimension of its character- 
istic variety equals n (i.e. dim(M) = n). An ideal / in Dx is called 
holonomic if either / = Dx or Dx/I is holonomic. 

The characteristic cycle of a holonomic -Dj<— module M is by defini- 
tion the sum 

CChar(M) := m c C 
C irreducible component of Char(M) 

viewed as a cycle in the cotangent bundle T*C n . Here mc is the mul- 
tiplicity of C in Char(M). 

Holonomic modules are particularly nice from the point of view of 
computations due to the following fact. 

Theorem 2 (Stafford [25]). Every holonomic Dx-module is cyclic. 

Every holonomic module M can be thought of as M = D x -Q, where 
Q is a cyclic generator, hence, M = D x / P^m Dx {Q). 

One example of a holonomic module is R = C[x]; it is generated 
by 1 G R. Another is its localization Rf = C[x, as stated by a 
theorem by J. Bernstein. 

Theorem 3 ([21 [3]). Let f G C[x] be a nonzero polynomial. Then the 
Dx-module Rf is holonomic. Moreover, for some negative integer a 
the element f a G Rf generates Rf as a Dx-module. 

Computing the annihilator Ann^ x (/ a ) is the main topic of this ar- 
ticle. 

By [21 E] the largest possible exponent a = a(f) above is equal to 
the smallest integer root of the Bernstein- Sato polynomial b /(s) G Q[s], 
defined as the monic nonzero polynomial b(s) of the smallest possible 
degree such that 

b(s)f s = Q(s, x, d) ■ f s+ \ where Q G D x [s] := C[s] ® c D x . 

The idea of many localization algorithms (first formulated by Oaku [19] 
in 1996) is simple to state: 

(1) Compute Ann(/ S ) as an ideal of 

(2) Compute bf(s) to determine a = a(f). 

(3) Specialization: Ann(/ S )| s=a . 

This approach, for example, is implemented in the function Dlocalize 
of the Dmodules package. The bottlenecks of the algorithm are items 
1 and 2 that require an expensive elimination via Grobner bases in the 
Weyl algebra Dx- 

Step 2 can be avoided by using the following estimate 
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Theorem 4 ([H], see also [21], [22] )• Rf is generated by f~ n+1 . 

As a corollary, for planar curves / e C[x, y], localization Rf = Dx • 
f-^Dx/Annif- 1 ). 

Remark 1. Note that one can easily compute the annihilator of a 
rational function in the the algebra of linear differential operators with 
rational function coefficients C(a;)[9] = <C(x)®c\x\Dx : the annihilator 
is a maximal ideal of C(x)[d], e.g., 

df_ 

dxi 



Ann cixm (f a ) = C(x) ® Ann Dx (f a ) = (fdt- a^\i = 1,. . . ,n). 



2. Iterative algorithm 

Our main example will be the family of the so-called Reiffen's curves; 
see [22] . 

f = f P , q = x p + y q + xy q - 1 = 0, 
where p > 4, q > p + 1. 

Example 1. Computing the annihilator (this can be done via AnnFs 
followed by specialization s = — 1), 

Ann(/ 4 ~ 5 1 ) = (Ax 2 d x + 5xyd x + 3xyd y + Ay 2 d y + 16x + 20y, 
16xy 2 d x + Ay 3 d x + 12y 3 d y - \2hxyd x - Ax 2 d y + 

hxyd y - 100y 2 d y + My 2 - 500y, 
16y 3 d 2 x - 16y 3 d x d y + 125xyd 2 x - 35xyd x d y + 100y 2 d x d y + 

12x 2 8 2 - 2xyd 2 - 24y 2 d 2 + \\2xyd x - m V 2 d x + 

My 2 d y - 930x9 x + 625^. + 2Qxd y - 

893j/0 tf + 448y-372O), 

we see that it is generated in order 2, i.e., generated by operators of 
order at most 2. 

One interesting question is: What is the order of generation of this 
annihilator for a given p and q ? 

2.1. Iterative approach. 

Definition 2. We call 

Ann (d) (/ a ) = (Qe Ann(/ a ) | ordQ < d) 

the d-th truncated annihilator of f a . 

To compute Ann^(/ a ) one may find the i?-syzygy module Sj for 
the vector of partial derivatives (d a ■ f a ) a where |a| < d. 
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Example 2. To find an element mAnn' 1 ^/ 1 ) for f = x 2 —y 3 consider 
all partial derivatives of f^ 1 of order at most 1: 

dx-r 1 = -2x r 2 

dyf- 1 = 3y 2 f~ 2 
1 • / ' = (r 2 - !!■'') / 2 

For instance, 3x(— 2x) + 2y(3y 3 ) + 6(x 2 — y 3 ) = 0, hence, 

3xd x + 2yd y + 6 G Ann (1) (/ -1 ). 

In Dmodules the function for computing the truncated annihilator is 
called kOrderAnnFa: 

11 : loadPackage "Dmodules"; 

12 : R = QQ[x,y] ; 

13 : f = x~2-y~3; 

14 : Al = kOrderAnnFa(l,f ,-1) 

2 3 2 2 

o4 = ideal (3x*dx + 2y*dy +6, 3y dx + 2x*dy, y dy - x dy + 3y ) 

Note that kOrderAnnFa(d,f ,-1) would return the same ideal for all 

d G N as the annihilator for this particular example is generated in 

order 1. 

Remark 2. Suppose S^-i, or in other words Aim.( d ~ 1 \f a ) , is known. 
Then the computation of Sj, or the d-th truncated annihilator, can be 
optimized by computing the syzygies modulo S^-i- 

In practice, however, it appears that the time of the syzygy computa- 
tion at order d dominates that for order d — 1 for all d. In particular, 
in Experiment Q] for the Reiffen curve of any degree the last step takes 
more time than all previous steps combined. 

2.2. Stopping criterion. The sequence of truncated annihilators sta- 
bilizes at some point: 

Ann (1) (/ a ) C Ann (2) (/ a ) C • • • C Ann (d) (/ a ) = Ann (d+1) (/ a ) 

Definition 3. The smallest d such that Ann^(/°) = Ann(/ a ) is de- 
noted by n(f a ) and is called the annihilator order of f a . 

Note that Ann (d) (/ a ) = Ann (d - 1) (/ a ) does not imply Ann (d+1) (/ a ) = 
Ann^(/ a ). For example, for / = x and a = 3 we have Ann^ = 
Anrv® = Ann^ = (xd x — 3) and Ann^ = (xd x — 3, d^). 

2.3. Annihilator order of a planar curve. In view of Theorem [4j 
for a planar curve defined by / G R := C[x, y] the localization Rf is 
generated by We have constructed an algorithm that computes 
truncated annihilators stopping exactly at the annihilator order of f~ l 
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thus producing the whole annihilator Ann(/ x ) for a plane curve with 
at most one singular point, that we can assume to be the origin. 

Algorithm 1. (d,A) =KappaAndAnnihilator(/) 

Input: f e R = C[x, y], a plane curve with at most one singular point, 

that we can assume to be the origin. 
Output: d = n(f~ l ), A = Ann(/^ 1 ). 

1: d<- 0. 

2: repeat 

3: d<-d+l. 

4: A 4- Ann^(/ _1 ). 

5: o primary component of gr w (A) corresponding to the origin. 

6: until dega = [multiplicity of f at the origin) — 1 

Proof of correctness. We see X = C 2 as a complex manifold and con- 
sider C, the analytic plane curve defined by the polynomial equation 
/ = 0. We assume / to be a reduced polynomial. Denote 

• by Ox the sheaf of holomorphic functions on X, 

• by T>x the sheaf of linear partial differential operators on X 
with holomorphic coefficients, and 

• by Ox [*C] the sheaf of meromorphic functions on X with poles 
on C. 

By [H] the sheaf Ox[*C] is a coherent holonomic P^-module. If p G 
X \ C then is just Ox, P and if the point p G C is a smooth 

point then the Px,p~ m odule Ox [*C] P is isomorphic to the quotient 

Pi,(o,o) 

^x,(o,o)(xd x + l,d y ) 

(for some choice of local coordinates (x,y)). So, in the neighborhood 
of a smooth point in C, the characteristic variety of is just 

the union of the conormal space to C, say TqX and the zero section 
T X X. Moreover, the multiplicities of each of these components in the 
corresponding characteristic cycle is 1. 

It is well known that (see e.g. [51 Section 6]) in a sufficiently small 
neighborhood U of a point p G C, the characteristic cycle of 0[/[*C] is 

CChar(O c/ [*C]) = T&J + T* nC ^ {p} U + (m - 1)T;U 

where TyU is the zero section of the cotangent bundle T*U, T UnC ,j p yU 
stands for the closure in T*U of the conormal bundle to the smooth 
part of U n C, T*U is the conormal bundle to the point p and m is the 
multiplicity of the plane curve C at the point p. 
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Let us denote X := Ann Vu (f 1 ) and for any integer number d > 1, 
:= Ann$ (f~ x ) considered as a sheaf of ideals in Vjj. 



As / is a global equation for the plane curve Cclwe have X 



(d) 
■q 

V u ^ q Ann Dx (f- 1 ) for any q G U. 

We have the following exact sequence of holonomic ©{/-modules: 

— >■ 77777 — 77777 — —^r- — >■ 

and then the following equality of characteristic cycles: 
bar I ^ 1 = CChar [ — 1 



Denote by the multiplicity of the irreducible component T*U 
in CChar (|py) • As the ideals X^ and X coincide on £/ fl C \ {p} the 
holonomic TV-module 7J3J * s concentrated on p. Then CChar (^%) = 
CChar (^-) if and only mS d ^ — m — 1 (where m is the multiplicity of 
the plane curve C at p E C). We also have X^^ = X if and only if 

m (d) _ m _ l 

Once we have a finite set of polynomial equations for the character- 
istic ideal gr w (Ann^ , these equations define the characteristic 
cycle CChar (|%) in the neighborhood U of any point p E X. The 
multiplicity equals the degree of the primary ideal in the primary 
decomposition of gr w (Ann^ corresponding to T*U. □ 

Remark 3. The multiplicity can be also obtained as the multiplic- 
ity of the (unique) pointp E T*U in the intersection gr w (Ann( d > (f~ 1 ))n 
L where L is the defining ideal of a generic 2-plane in the 4-dimensional 
ambient space. 

Remark [3] allows us to avoid the expensive primary decomposition 
step in the Algorithm[TJ Assume p = is the only singularity; hence, let 
U = X and let S = C[x, y, £, 77] be the ring of the characteristic ideal 
gr w {Ann^ d \f^ 1 )). In our implementation we set L = (£ — a, n — b), 
(a, b) E C 2 \ {(0, 0)}; in practice, a and b are taken to be small integers. 
The genericity of the plane L means that the point p = (0, 0, a, b) does 
not belong to X^ nC « f p y U. The following algorithm checks the genericity 
of the choice. 

Algorithm 2. g =CheckGenericity(/, a, b) 

Input: f G R = C[x,y], a plane curve with the unique singular point 
at the origin. 

Output: g, a boolean value, whether the genericity condition described 
in Remark^ is satisfied. 
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1: Compute the syzygies M <— {(u, v) G R 2 \ uf x + vf y = 0}. 

2: I — {u£ + VT) I («, v) G 5) + (/) . 

5: Compute the saturation J <— 7 : (x,y)°°. 

4: g 4— (0, 0, a, 6) £ V( J) 

Remark 4. TVie choice of the plane L = (£ — a, 77 — 6), (a, 6) G C 2 \ 
{(0,0)}, is artificial but sufficient, since dim(T*U PI T^ nCr > r > U) < 1. 
In /act, in Experiment^ the choice (a, 6) = (0, 1) always works. 

Remark 5. In view of Remark^ can be computed numerically: 

one can pick a random L and compute an approximation to the point 
of intersection p. Using this approximation, one may recover the mul- 
tiplicity heuristically . To the best of our knowledge there is no 
known approach to certification of the correctness of such a computa- 
tion. 

Moreover, in our experiments the bottleneck is not the saturation 
step in Algorithm^ but the computation of syzygies (described after 
Definition^). 

Remark 6. The algorithm can be generalized for arbitrary complex 
planar curves C defined by a reduced polynomial f . As the number of 
singular points ofC is finite it is enough to apply our iterative algorithm 
at each singular point of C and compute ^{fp 1 ) for each p G Sing(C). 
We finally have n(f~ l ) as the maximum of the ^(fp 1 ) forp G SingiC). 
See also Section \3.1\ for a generalization of this algorithm to a more 
general situation. 

We have implemented Algorithm [1] modified according to the Re- 
mark [3] and including the genericity check provided by Algorithm [2] in 
the Dmodules package for Macaulay2. The following example demon- 
strates that the new function kappaAnnFlPlanarCurve produces the 
same result as the old method involving AnnFs. 

12 : f = reiff en(7,8) 

7 8 7 
o2=xx +x +x 
12 2 1 

02 : QQ[x , x ] 

1 2 

13 : As = AnnFs f ; 

03 : Ideal of QQ [x , x , dx , dx , s] 

12 12 

14 : A = sub(As, {last gens ring As => -1}); 

15 : (kappa, A') = kappaAnnFlPlanarCurve f ; 

16 : kappa 
06 : 4 
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i7 : A == sub (A' , ring A) 
o7 = true 

Experiment 1. We have computed the annihilator order for the Reif- 
fen curve: K^f~V) seems to depend only on p. Starting with p = 4 (the 
last computed value is forp = 21), the sequence n(f~l) is 

2, 2, 3, 4, 4, 5, 6, 6, 7, 8, 8, 9, 10, 10, 11, 12, 12, 13, . . . 

We can compute K>(fpq) for q < 22 in less than one day using our 
implementation of Algorithm 1. In contrast, an implementation of the 
general algorithm (the computation o/Ann(/ s ) followed by the special- 
ization s = — 1) took about a week for p = 17 and has not produced a 
result in one month forp = 18. 



Another interesting observation is that the sequence of multiplicities 
(m^) for d — 1, . . . , K(fpq) does not depend on q. Below we report 
the computed values rrv®: 
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3. Discussion and open problems 

We have proposed a new algorithm computing a presentation, as 
a left A 2 -module, of the localization C[x, y]f for nonzero polynomials 
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/ = f(x,y) in two variables. The proof of the correctness of the algo- 
rithm is specific to dimension 2 and it seems difficult to generalize it to 
higher dimensions. Nevertheless, in the next subsection we describe an 
analogous algorithm for polynomials / = f(xi, . . . , x n ), for any n > 0, 
under the hypothesis that the hypersurface defined by / in C n has only 
isolated singularities. 

3.1. Isolated hypersurface singularities. The iterative algorithm 
described before can be also applied in a more general situation. Let 
n > be an integer. Denote X = C n , by Ox the sheaf of rigs 
of holomorphic functions on X and by T>x the sheaf of rings of lin- 
ear differential operators with holomorphic coefficients. Assume that 
/ G R = C[x±, . . . ,x n ] is a non zero reduced polynomial defining a 
hypersurface C G X with an isolated singularity p G C. Then the 
characteristic cycle of 0x[*C] has the following description (see e.g. 
[5j Section 6]) 



CChax(O x [*C]) = T X X + T* X{p] X + vT;X 

where v is the Milnor number of a generic hyperplane section of C 
through the point p. This formula coincides with the one we have 
considered before for n = 2. The number v can be computed in a 
effective way. We can repeat the same ideas as in the proof of the 
correctness of Algorithm [TJ By comparison of the multiplicity rrv& 
with v, where is the multiplicity of the irreducible component 

t;x in 



CChar 



V x Ann An (f- n+1 
we can compute the annihilator order of f~ n+1 . 

3.2. Weyl closure. 

Definition 4. The Weyl closure of I C Dx is defined as 
WeylCl(J) = K{x) ® c[x] lnD x . 

Note that the problem of computing Ann(/ a ) is a particular instance 
of the problem of computing the Weyl closure: according to [1] the 
annihilator of f a in K(x) ®c{x] Dx is equal to (fdi + J^- | i — 1, . . . , n). 
Therefore, 

Ann(r 1 ) = WeylCl«/9 i + J£ | % = 1, . . . , n)). 

There exists an algorithm due to Tsai [28] that computes the Weyl 
closure. However, the computation is highly nontrivial and boils down 
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to computing localizations of Dx-modules anyway. On the other hand, 
one may define the truncated Weyl closure of order d analogously to the 
d-th truncated annihilator. Then, the sequence of truncated closures 
stabilizes at some degree; finding this degree or discovering a termina- 
tion criterion for an iterative algorithm similar to Algorithm [1] would 
result in a much simpler procedure. 
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